########################################################################
# PLS 206 — Week 3 live coding, part 1: groups as predictors
# Grey Monroe, Fall 2026
########################################################################

# install.packages(c("ggplot2", "emmeans", "car"))
library(ggplot2)

########################################################################
# 1. The data
########################################################################

# 90 plots: 3 varieties x 2 water treatments, nitrogen varies by plot.

trial <- read.csv("https://greymonroe.github.io/PLS_206/data/wheat_trial.csv")
head(trial)
table(trial$variety, trial$water)

ggplot(trial, aes(water, yield)) +
  geom_boxplot() +
  geom_jitter(width = 0.1) +
  theme_classic()

########################################################################
# 2. Two groups: a dummy variable
########################################################################

# Intercept = mean of the reference group (control).
# waterdrought = drought minus control.

m_water <- lm(yield ~ water, data = trial)
summary(m_water)
tapply(trial$yield, trial$water, mean)

# the same test as a t-test
t.test(yield ~ water, data = trial, var.equal = TRUE)

# the 0/1 column R makes for you
head(model.matrix(m_water))

# choose the reference level
trial2 <- trial
trial2$water <- factor(trial2$water, levels = c("drought", "control"))
coef(lm(yield ~ water, data = trial2))

########################################################################
# 3. Three groups: ANOVA
########################################################################

control <- subset(trial, water == "control")
m_var <- lm(yield ~ variety, data = control)
summary(m_var)      # each row: that variety minus Amber
anova(m_var)        # do the varieties differ at all?
TukeyHSD(aov(m_var))   # which pairs?

# effect size: share of variation explained by variety
ss <- anova(m_var)$"Sum Sq"
ss[1] / sum(ss)

########################################################################
# 4. Two factors and their interaction
########################################################################

ggplot(trial, aes(variety, yield, colour = water)) +
  geom_boxplot() +
  theme_classic()

m_add <- lm(yield ~ variety + water, data = trial)   # one drought effect
m_int <- lm(yield ~ variety * water, data = trial)   # one per variety
anova(m_int)        # test the interaction row first
summary(m_int)

# the drought effect in each variety
library(emmeans)
pairs(emmeans(m_int, ~ water | variety))

########################################################################
# 5. summary() vs anova() with an interaction
########################################################################

# varietyCopper = Copper vs Amber in the reference water level only.
# Change the reference: the coefficient changes, anova() doesn't.

m_int2 <- lm(yield ~ variety * water, data = trial2)
coef(summary(m_int))["varietyCopper", ]    # under control
coef(summary(m_int2))["varietyCopper", ]   # under drought
anova(m_int)
anova(m_int2)

########################################################################
# 6. Unbalanced data: ANOVA types
########################################################################

unbal <- trial[-(1:20), ]   # lose 20 plots
anova(lm(yield ~ variety + water, data = unbal))
anova(lm(yield ~ water + variety, data = unbal))   # order matters

library(car)
Anova(lm(yield ~ variety + water, data = unbal), type = 2)

########################################################################
# 7. A group and a slope: ANCOVA
########################################################################

m_anc <- lm(yield ~ nitrogen + water, data = trial)    # parallel lines
summary(m_anc)
m_anc2 <- lm(yield ~ nitrogen * water, data = trial)   # different slopes
summary(m_anc2)

ggplot(trial, aes(nitrogen, yield, colour = water)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic()

########################################################################
# 8. Two numeric predictors: centre before reading main effects
########################################################################

# nitrogen and water (mm/week) both numeric (simulated). With
# nitrogen * water, the nitrogen row is the nitrogen slope at water = 0.

set.seed(47)
irrig <- data.frame(nitrogen = round(runif(60, 0, 200)),
                    water = round(runif(60, 15, 45)))
irrig$yield <- round(2 + 0.04 * irrig$water +
  0.0004 * irrig$nitrogen * irrig$water + rnorm(60, 0, 0.6), 2)

m_nw <- lm(yield ~ nitrogen * water, data = irrig)
summary(m_nw)        # nitrogen: nothing?
anova(m_nw)          # nitrogen: strong

irrig$nitrogen_c <- irrig$nitrogen - mean(irrig$nitrogen)
irrig$water_c    <- irrig$water - mean(irrig$water)
m_nwc <- lm(yield ~ nitrogen_c * water_c, data = irrig)
summary(m_nwc)       # nitrogen slope at average water